module ghs
implicit none
contains
function ghcfun0(beta,eta,si)
use bessels
double precision beta(:), eta,si,nu,gam
double precision ghcfun0,a,raiz,x

nu=-1/(2*eta)
gam=(1-si)/si


a=besseld(nu+1.0d0,gam)-1.0d0
x=dot_product(beta,beta)
if (x >1.0d-4) then
    raiz=dsqrt(1.0d0+4.0d0*a*x)
    ghcfun0=(-1.0d0+raiz)/(2.0d0*a*x)
else
    ghcfun0=1.0d0-a*x+2.0d0*((a*x)**2.0d0)-5.0d0*((a*x)**3.0d0)+14.0d0*((a*x)**4.0d0)
endif


end function ghcfun0

function ghcfun1(bb,eta,si)
use bessels
double precision bb, eta,si,nu,gam
double precision ghcfun1,a,raiz,x

nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si


a=besseld(nu+1.0d0,gam)-1.0d0
x=bb
if (x >1.0d-4) then
    raiz=dsqrt(1.0d0+4.0d0*a*x)
    ghcfun1=(-1.0d0+raiz)/(2.0d0*a*x)
else
    ghcfun1=1.0d0-a*x+2.0d0*((a*x)**2.0d0)-5.0d0*((a*x)**3.0d0)+14.0d0*((a*x)**4.0d0)
endif


end function ghcfun1


subroutine ghparams0(alpha,upsilon,beta,eta,si)
use bessels
use linear_operators
double precision, intent(out):: alpha(:), upsilon(:,:)
double precision, intent(in):: beta(:),eta,si
double precision mbeta(1:size(beta),1:1),nu,gam,coef,c
integer n

n=size(beta)
c=ghcfun0(beta,eta,si)


coef=((c-1)/dot_product(beta,beta))
mbeta=matvec(beta)

nu=-1/(2*eta)
gam=(1-si)/si

alpha=-c*beta
upsilon=(eye(n)+coef*matmul(mbeta,transpose(mbeta)))*gam/besselr(nu,gam)
end subroutine ghparams0

function stdghdensity(y,beta,eta,si,coef,besrat,iomegtstar,alpha,cfun)
!GH density function  of a standardised GH vector,constant term comes from outside
use bessels
use linear_operators
double precision y(:),beta(:),eta,si
double precision nu,gam,lt(1:size(y)),arg
double precision stdghdensity,besrat,iomegtstar(1:size(y),1:size(y)),qfun,coef,a&
                 &,alpha,cfun
integer n

n=size(y)
nu=-1.0d0/(2.0d0*eta)
gam=(1-si)/si

lt=y+cfun*beta
qfun=dsqrt(1.0d0+besrat*(dot_product(lt,matmul(iomegtstar,lt))))


arg=alpha*qfun

stdghdensity=dexp(coef+(nu-.5*n)*dlog(qfun)+robustlbesk(nu-.5*n,arg)+dot_product(beta,lt))

end function stdghdensity
!!!!!!!!!!
subroutine stdghcons(coef,besrat,iomegtstar,alpha,cfun,beta,eta,si)
!Constant term of GH log-density function
use bessels
use linear_operators
double precision, intent(in)::beta(:),eta,si
double precision, intent(out)::coef,besrat,iomegtstar(:,:),alpha,cfun
double precision nu,gam,normb,mbeta(1:size(beta),1:1),arg
double precision omegtstar(1:size(beta),1:size(beta)),a
integer n


n=size(beta)
nu=-1.0d0/(2.0d0*eta)
gam=(1-si)/si

cfun=ghcfun0(beta,eta,si)
besrat=besselr(nu,gam)/gam
normb=dot_product(beta,beta)
alpha=dsqrt((1/besrat)*cfun*normb+gam**2.0d0)
mbeta=matvec(beta)
if (normb .gt. 1.0d-3) then
	iomegtstar=eye(n)-((cfun-1.0d0)/(cfun*normb))*(mbeta .xt. mbeta)
	omegtstar=eye(n)+((cfun-1.0d0)/normb)*(mbeta .xt. mbeta)
else
	a=besseld(nu+1.0d0,gam)-1.0d0
	iomegtstar=eye(n)-((-a+2.0d0*(a*normb)**2.0d0-5.0d0*a**3.0d0*normb**2.0d0)/cfun)*(mbeta .xt. mbeta)
	omegtstar=eye(n)+(-a+2.0d0*(a*normb)**2.0d0-5.0d0*a**3.0d0*normb**2.0d0)*(mbeta .xt. mbeta)
endif

coef=nu*dlog(gam)-.5*n*dlog(2.0d0*3.14159265358979d0)-.5d0*(nu-.5d0*n)*dlog(gam*gam+(cfun*normb/besrat))&
    &-.5d0*dlog(det(omegtstar))+.5d0*n*dlog(besrat)-robustlbesk(nu,gam)
end subroutine stdghcons



!!!!!!!!!!
function matvec(x)
double precision x(:)
integer n,i
double precision matvec(1:size(x),1:1)

n=size(x)

do i=1,n,1
	matvec(i,1)=x(i)
end do
end function matvec
!!!!!!!!!!!!!
subroutine ghskku(sk,ku,beta,eta,si)
use bessels
double precision beta,vbeta(1), eta,si,sk,ku,coef3,coef4,nu,gam,cfun,d,c
nu=-1.0d0/(2.0d0*eta)
gam=(1.0d0-si)/si

coef3=besselr(nu+2.0d0,gam)*besselr(nu+1.0d0,gam)/(besselr(nu,gam)**2.0d0)
coef4=besselr(nu+3.0d0,gam)*besselr(nu+2.0d0,gam)*besselr(nu+1.0d0,gam)/(besselr(nu,gam)**3.0d0)
vbeta(1)=beta
cfun=ghcfun0(vbeta,eta,si)
d=besseld(nu+1.0d0,gam)

sk=cfun**3.0d0*beta**3.0d0*(coef3-3*d+2.0d0)&
&+3*cfun**2.0d0*beta*(d-1)
ku=cfun**4.0d0*beta**4.0d0*(coef4-4.0d0*coef3+6.0d0*d-3.0d0)&
&+6*cfun**3.0d0*beta**2.0d0*(coef3-2*d+1.0d0)+3*cfun**2.0d0*d
end subroutine ghskku

subroutine marginals(beta1,w,var,b,eta,si)
use linear_operators
double precision beta1,beta2,var(1:2,1:2),b(2),eta,si,&
&nu,gam,beta(1:2),cfun,romt(1:2,1:2),w(1:2)

romt=.t.chol(var)
beta=(.t.romt).x.b
cfun=ghcfun0(beta,eta,si)
if (dot_product(b,b)==0.0d0) then
	beta1=0.0d0
else
	beta1=cfun*dot_product(w,beta)*dsqrt(dot_product(w,w))/(dot_product(w,w)+(cfun-1.0d0)*(dot_product(w,beta)**2.0d0/dot_product(beta,beta)))
endif
end subroutine marginals
end module ghs
